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Self-Organized Criticality is the emergence of long-ranged spatio-temporal correlations in non- 
equilibrium steady states of slowly driven systems without fine tuning of any control parameter. 
Sandpiles were proposed as prototypical examples of self-organized criticality. However, only some 
of the laboratory experiments looking for the evidence of criticality in sandpiles have reported a 
positive outcome. On the other hand a large number of theoretical models have been constructed 
that do show the existence of such a critical state. We discuss here some of the theoretical models 
as well as some experiments. 



The concept of Self-Organized Criticality (SOC) was 
introduced by Bak, Tang and Wiesenfeld (BTW) in 1987 
0. It says that there is a certain class of systems in 
nature whose members become critical under their own 
dynamical evolutions. An external agency drives the sys- 
tem by injecting some mass (in other examples, it could 
be the slope, energy or even local voids) into it. This 
starts a transport process within the system: Whenever 
the mass at some local region becomes too large, it is 
distributed to the neighbourhood by using some local re- 
laxation rules. Globally, mass is transported by many 
such successive local relaxation events. In the language 
of sandpiles, these together constitute a burst of activity 
called an avalanche. If we start with an initial uncritical 
state, initially most of the avalanches are small, but the 
range of sizes of avalanches grows with time. After a long 
time, the system arrives at a critical state, in which the 
avalanches extend over all length and time scales. Cus- 
tomarily, critical states have measure zero in the phase 
space. However, with self-organizing dynamics, the sys- 
tem finds these states in polynomial times, irrespective 
of the initial state (2|-||] . 

BTW used the example of a sandpile to illustrate their 
ideas about SOC. If a sandpile is formed on a horizon- 
tal circular base with any arbitrary initial distribution 
of sand grains, a sandpile of fixed conical shape (steady 
state) is formed by slowly adding sand grains one after 
another (external drive). The surface of the sandpile in 
the steady state on the average makes a constant angle 
known as the angle of repose, with the horizontal plane. 
Addition of each sand grain results in some activity on 
the surface of the pile: an avalanche of sand mass fol- 
lows, which propagates on the surface of the sandpile. 
Avalanches are of many different sizes and BTW argued 
that they would have a power law distribution in the 
steady state. 

There are also some other naturally occurring phe- 
nomena which are considered to be examples of SOC. 
Slow creeping of tectonic plates against each other re- 
sults intermittent burst of stress release during earth- 



quakes. The energy released is known to follow power law 
distributions as described by the well known Gutenberg- 
Richter Law The phenomenon of earthquakes is be- 
ing studied using SOC models Q. River networks have 
been found to have fractal properties. Water flow causes 
erosion in river beds, which in turn changes the flow dis- 
tribution in the network. It has been argued that the 
evolution of river pattern is a self-organized dynamical 
process 0. Propagation of forest fires j|] and biologi- 
cal evolution processes |^ have also been suggested to be 
examples of SOC. 
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FIG. 1. An avalanche of the Abelian Sandpile Model, 
generated on a 3 x 3 square lattice. A sand grain is dropped 
on a stable configuration at the central site. The avalanche 
created has size s = 6, area a = 6, life-time t = A and the 
radius r = \f2. 

Laboratory experiments on sandpiles, however, have 
not always found evidence of criticality in sandpiles. In 
the first experiment, the granular material was kept in 
a semicircular drum which was slowly rotated about the 
horizontal axis, thus slowly tilting the free surface of the 
pile. Grains fell vertically downward and were allowed to 
pass through the plates of a capacitor. Power spectrum 
analysis of the time series for the fluctuating capacitance 
however showed a broad peak, contrary to the expecta- 
tion of a power law decay, from the SOC theory jio). 

In a second experiment, sand was slowly dropped on 
to a horizontal circular disc, to form a conical pile in the 
steady state. On further addition of sand, avalanches 
were created on the surface of the pile, and the outflow 
statistics was recorded. The size of the avalanche was 
measured by the amount of sand mass that dropped out 
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of the system. It was observed that the avalanche size 
distribution obeys a scahng behaviour for small piles. For 
large piles, however, scaling did not work very well. It 
was suggested that SOC behavior is seen only for small 
sizes, and very large systems would not show SOC Q. 

Another experiment used a pile of rice between two ver- 
tical glass plates separated by a small gap. Rice grains 
were slowly dropped on to the pile. Due to the anisotropy 
of grains, various packing configurations were observed. 
In the steady state, avalanches of moving rice grains re- 
freshed the surface repeatedly. SOC behaviour was ob- 
served for grains of large aspect ratio, but not for the less 
elongated grains |p^ . 

Theoretically, however, a large number of models have 
been proposed and studied. Most of these models study 
the system using cellular automata where discrete, as 
well as continuous, variables are used for the heights of 
sand columns. Among them, the Abelian Sandpile model 
is most popular ||l|,|l^. Other models of self organized 
criticality have been studied but will not be discussed 
here. These include the Zheng model which has modi- 
fied rules for sandpile evolution |0 , a model for Abelian 
distributed processors and other stochastic rule models 
the Eulerian Walkers model and the Takayasu 
aggregation model . 

In the Abelian sandpile model, we associate a non- 
negative integer variable h representing the height of the 
'sand column' with every lattice site on a d-dimensional 
lattice (in general on any connected graph). One often 
starts with an arbitrary initial distribution of heights. 
Grains are added one at a time at randomly selected sites 
O: ho — > ho + 1. The sand column at any arbitrary site 
i becomes unstable when hi exceeds a previously selected 
threshold value he for the stability. Without loss of gen- 
erality, one usually chooses he = 2d — 1. An unstable 
sand column always topples. In a toppling, the height 
is reduced as: hi ^ hi — 2d and all the 2d neighbouring 
sites {j} gain a unit sand grain each: hj —>■ hj + 1. This 
toppling may make some of the neighbouring sites unsta- 
ble. Consequently, these sites will topple again, possibly 
making further neighbours unstable. In this way a cas- 
cade of topplings propagates, which finally terminates 
when all sites in the system become stable (Fig. 1). One 
waits until this avalanche stops before adding the next 
grain. This is equivalent to assuming that the rate of 
adding sand is much slower than the natural rate of re- 
laxation of the system. The wide separation of the 'time 
scale of drive' and 'time scale of relaxation' is common in 
many models of SOC. For instance, in earthquakes, the 
drive is the slow tectonic movement of continental plates, 
which occurs over a timescale of centuries, while the ac- 
tual stress relaxation occurs in quakes, whose duration is 
only a few seconds. This separation of time scales is usu- 
ally considered to be a defining characteristic of SOC. 
However, Dhar has argued that the wide separation of 
time scales should not be considered as a necessary con- 



dition for SOC in general Q. Finally, the system must 
have an outlet, through which the grains go out of the 
system, which is absolutely necessary to attain a steady 
state. Most popularly, the outlet is chosen as the {d — 1) 
dimensional surface of a d-dimensional hypercubic sys- 
tem. 
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FIG. 2. An example to show that a directed slope model 
is non- Abelian. Two slopes are measured from any site (i, j) 
as h{i,j) — h{i,j + 1) and h{i,j) — h{i -\- l,j -\- 1). If either of 
them is greater than 1, two grains are transferred from (i,j) 
and are given one each to ~\- 1) and (i -|- l,j -|- 1). On 
dropping a grain on the initial stable configuration, we see 
that finally two difi'erent height configurations result due to 
two difi'erent sequences of topplings [20]. 

The beauty of the Abelian model is that the final sta- 
ble height configuration of the system is independent of 
the sequence in which sand grains are added to the sys- 
tem to reach this stable configuration (isj. On a stable 
configuration C, if two grains are added, first at i and 
then at j, the resulting stable configuration C is exactly 
same in case the grains were added first at j and then 
at i. In other sandpile models, where the stability of 
a sand column depends on the local slope or the local 
Laplacian, the dynamics is not Abelian, since toppling of 
one unstable site may convert another unstable site to a 
stable site (Fig. Many such rules have been studied 
in the literature [l^Jl^ . 

An avalanche is a cascade of topplings of a number 
of sites created on the addition of a sand grain. The 
strength of an avalanche in general, is a measure of the 
effect of the external perturbation created due to the ad- 
dition of the sand grain. Quantitatively, the strength of 
an avalanche is estimated in four different ways: (i) size 
(s): the total number topplings in the avalanche, (ii) area 
(a): the number of distinct sites which toppled, (iii) life- 
time {t): the duration of the avalanche and (iv) radius 
(r) : the maximum distance of a toppled site from the ori- 
gin. These four different quantities arc not independent 
and are related to each other by scaling laws. Between 
any two measures x,y G {s,o,i,r} one can define a mu- 
tual dependence as: < y >~ x^^". These exponents are 
related to one another, e.g., ^ts = Itrlrs- For the ASM, 
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it can be proved that the avalanche clusters cannot have 
any holes. It has been shown that ^rs = 2 in two dimen- 
sions. It has also benn proved that 7^4 = 5/4 A 
better way to estimate the ^tx exponents is to average 
over the intermediate values of the size, area and radius 
at every intermediate time step during the growth of the 
avalanche. 

Quite generally, the finite size scaling form for the 
probability distribution function for any measure x £ 
{s, a, i, r} is taken to be: 



P{x) 



The exponent Ux determines the variation of the cut-off 
of the quantity x with the system size L. Alternatively, 
sometimes it is helpful to consider the cumulative prob- 
ability distribution F{x) — P{x)dx which varies as 
^.i-T^^ However, in the case of Tx = 1, the variation 
should be in the form F{x) = C— log(a;). Between any 
two measures, scaling relations like jxy = ('''x — 1)/(tj, — 1) 
exist. Recently, the scaling assumptions for the avalanche 
sizes have been questioned. It has been argued that there 
actually exists a multifractal distribution instead |2l| . 

Numerical estimation for the exponents have yielded 
scattered values. For example estimates of the exponent 
Ts range from 1.20 to 1.27 (H and 1.29 [||. 

We will now look into the structure of avalanches in 
more detail. A site i can topple more than once in the 
same avalanche. The set of its neighbouring sites {j}, 
can be divided into two subsets. Except at the origin 
O, where a grain is added from the outside, for a top- 
pling, the site i must receive some grains from some of 
the neighbouring sites {ji} to exceed the threshold he- 
These sites must have toppled before the site i. When 
the site i topples, it loses 2c? grains to the neighbours, by 
giving back the grains it has received from {ji}, and also 
donating grains to the other neighbours {^2}- Some of 
these neighbours may topple later, which returns grains 
to the site i and its height hi is raised. The following 
possibilities may arise: (i) some sites of {^2} may not 
topple at all; then the site i will never re-topple and is a 
singly toppled site on the surface of the avalanche, (ii) 
all sites in {^2} topple, but no site in {ji} topples again; 
then i will be a singly toppled site, surrounded by singly 
toppled sites, (iii) all sites in {^2} topple, and some sites 
of {ji} re-topple; then i will remain a singly toppled site, 
adjacent to the doubly toppled sites, (iv) all sites in {J2} 
topple, and all sites of {ji} re-topples; then the site i 
must be a doubly toppled site. This implies that the set 
of at least doubly toppled sites must be surrounded by 
the set of singly toppled sites. Arguing in a similar way 
will reveal that sites which toppled at least n times, must 
be a subset and also are surrounded by the set of sites 
which toppled at least (n — l) times. Finally, there will be 
a central region in the avalanche, where all sites have top- 
pled a maximum of m times. The origin of the avalanche 



O, where the sand grain was dropped, must be a site in 
this maximum toppled zone. Also the origin must be at 
the boundary of this to*^ zone, since otherwise it should 
have toppled (m + 1) times |2J]. 
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FIG. 8.- i-iaf Aahkas; nple -of 4heJlMiglit -distribution in a 
recurrent configuration C' on a 24 x 24 square lattice. This 
configuration is obtained by dropping a grain a some previous 
configuration C at the encircled site, (b) The spanning tree 
representation of the configuration C' . (c) A new configura- 
tion C" is obtained by taking out one grain at the encircled 
site from the configuration C' . A spanning tree cannot be ob- 
tained for C" . The bonds of the spanning tree corresonding to 
the forbidden sub-configuration in C" are shown by the thin 
lines. 

Using this idea, we see that the boundary sites on any 
arbitrary system can topple at most once in any arbitrary 
number of avalanches. Similar restrictions are true for 
inner sites also. A (2n-f- 1) x {2n+l) square lattice can be 
divided into (n + 1) subsets which are concentric squares. 
Sites on the m-th such square from the boundary can 
topple at most m times, where as the central site cannot 
topple more than n times in any avalanche. 

Avalanches can also be decomposed in a different way, 
using Waves of Toppling. Suppose, on a stable configu- 
ration C, a sand grain is added at the site O. The site is 
toppled once, but is not allowed to topple for the second 
time, till all other sites become stable. This is called the 
first wave. It may happen that after the first wave, the 
site O is stable; in that case the avalanche has termi- 
nated. If the site O is still unstable it is toppled for the 
second time, and all other sites are allowed to become 
stable again; this is called the second wave, and so on. It 
was shown, that in a sample where all waves occur with 
equal weights, the probability of occurrence of a wave of 
area a is D{a) ^ 1/a p5| ]. 

It is known that the stable height configurations in 



ASM are of two types: Recurrent configurations ap- 
pear only in the steady state with uniform probabili- 
ties, whereas Transzent configurations occur in the steady 
state with zero probability. Since long range correlations 
appear only in the steady states, it implies that the re- 
current configurations are correlated. This correlation 
is manifested by the fact that certain clusters of con- 
nected sites with some specific distributions of heights 
never appear in any recurrent configuration. Such clus- 
ters are called the forbidden sub-configurations. It is easy 
to show that two zero heights at the neighbouring sites: 
(0—0) or, an unit height with two zero heights at its two 
sides: (0—1—0) never occur in the steady state. There 
are also many more forbidden sub-configurations of big- 
ger sizes. 

An L X L lattice is a graph, which has all the sites 
and all the nearest neighbour edges (bonds). A Span- 
ning Tree is a sub-graph of such a graph, having all sites 
and some bonds. It has no loop and therefore, between 
any pair of sites there exists an unique path through a 
sequence of bonds. There can be many possible Spanning 
trees on a lattice. These trees have interesting statistics 
in a sample where they are equally likely. Suppose when 
we randomly select such a tree and then randomly se- 
lect one of the unoccupied bonds and occupy it, it forms 
a loop of length £. It has been shown that these loops 
have the length distribution D{() ^ Similarly, 
if a bond of a Spanning tree is randomly selected and 
deleted, then it divides into two fragments. The sizes of 
the two fragments generated follow a probability distri- 
bution D{a) ~ a~^^/* |26j. It was also shown that every 
recurrent configuration of the Abelian model on an arbi- 
trary lattice has a one-to-one correspondence to a random 
Spanning tree graph on the same lattice. Therefore, there 
are exactly the same number of distinct Spanning trees 
as the number of recurrent Abelian sandpile model con- 
figurations on any arbitrary lattice |^^. Given a stable 
height configuration, there exists an unique prescription 
to obtain the equivalent Spanning tree. This is called 
the Burning method [^o| . A fire front, initially at every 
site outside the boundary, gradually penetrates (burns) 
into the system using a deterministic rule. The paths of 
the fire front constitute the Spanning tree. A fully burnt 
system is recurrent, otherwise it is transient. 

Suppose, addition of a grain at the site O of a stable 
recurrent configuration C, leads to another stable config- 
uration C . Is it possible to get back the configuration 
C knowing C and the position of 07 This is done by 
Inverse toppling pTf . Since C is recurrent, a correspond- 
ing Spanning tree (C) exists. Now, one grain at O is 
taken out from C and the configuration C"= C — Soj is 
obtained. This means on ST(C'), one bond is deleted at 
O and it is divided into two fragments. Therefore one 
cannot burn the configuration C" completely since the 
resulting tree has a hole consisting of at least the sites of 
the smaller fragment (Fig. 3). This implies that C" has 
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a forbidden sub-configuration (Fi) of equal size and C" 
is not recurrent. On (Fi), one runs the inverse toppling 
process: 4 grains are added to each site i, and one grain 
each is taken out from all its neighbours {j}. The clus- 
ter of /i sites in Fi is called the first inverse avalanche. 
The lattice is burnt again. If it still has a forbidden 
sub-configuration (F2), another inverse toppling process 
is executed, and is called the second inverse avalanche. 
The size of the avalanche is: s = /i + /2 + /a + and 
the /i is related to the maximum toppled zone of the 
avalanche. From the statistics of random spanning trees 
p6[ it is clear that fi should have the same statistics of 
the two fragments of the tree generated on deleting one 
bond. Therefore the maximum toppled zone also has a 
power law distribution of the size, D{a) ~ a~^^/^. 

Sandpile models with stochastic evolution rules have 
also been studied. The simplest of these is a Two-state 
sandpile model, A stable configuration of this system con- 
sists of sites, either vacant or occupied by at most one 
grain. If there are two or more grains at a site at the same 
time we say there is a collision. In this case, all grains 
at that site are moved. Each grain chooses a randomly 
selected site from the neighbours and is moved to that 
site. The avalanche size is the total number of collisions 
in an avalanche. From the numerical simulations, the 
distribution of avalanche sizes is found to follow a power 
law, characterized by an exponent Ts ~ 1.27 This 
two-state model has a nontrivial dynamics even in one 
dimension . Recently, it has been shown that instead 
of moving all grains, if only two grains are moved ran- 
domly leaving others at the site, the dynamics is Abelian 

Some other stochastic models also have nontrivial crit- 
ical behaviour in one dimension. To model the dynam- 
ics of rice piles, Christensen et. al. studied the follow- 
ing slope model On a one-dimensional lattice of 
length L, non-negative integer variable hi represents the 
height of the sand column at the site i. The local slope 
Zi = hi ~ ft-i+i is defined, maintaining zero height on the 
right boundary. Grains are added only at the left bound- 
ary i = 1. Addition of one grain hi ^ hi + 1 implies an 
increase in the slope Zi ^ Zi + 1. If at any site, the lo- 
cal slope exceeeds a pre-assigned threshold value , one 
grain is transferred from the column at i to the column 
at (i + 1). This implies a change in the local slope as: 
Zi ^ Zi — 2 and Zi±i — > Zi±i -\- 1. The thresholds of the 
instability z^ are dynamical variables and are randomly 
chosen between 1 and 2 in each toppling. Numerically, 
the avalanche sizes are found to follow a power law dis- 
tribution with an exponent Ts « 1.55 and the cutoff ex- 
ponent was found to be cr^ « 2.25. This model is referred 
as the Oslo model. 

Addition of one grain at a time, and allowing the sys- 
tem to relax to its stable state, implies a zero rate of 
driving of the system. What happens when the driv- 
ing rate is finite? Corral and Paczuski studied the Oslo 



model in the situation of nonzero flow rate. Grains were 
added with a rate r, i.e., at every (1/r) time updates, 
one grain is dropped at the left boundary « = 1. They 
observed a dynamical transition separating intermittent 
and continuous flows [ p2[ . 

Many different versions of the sandpile model have 
been studied. However the precise classiflcation of var- 
ious models in different universality classes in terms of 
their critical exponents is not yet available and still at- 
tracts much attention [|l^,0 . Exact values of the criti- 
cal exponents of the most widely studied Abelian model 
are still not known in two dimensions. Some effort has 
also been made towards the analytical calculation of 
avalanche size exponents ^3|-^. Numerical studies for 
these exponents are found to give scattered values. On 
the other hand the two-state sandpile model is believed 
to be better behaved and there is good agreement of nu- 
merical values of its exponents by different investigators. 
However, whether the Abelian model and the two-state 
model belong to the same universality class or not is still 
an unsettled question 

If a real sandpile is to be modeled in terms of any of 
these sandpile models or their modifications, it must be 
a slope model, rather than a height model. However, 
not much work has been done to study the slope models 
of sandpiles (l^Jl^ . Another old question is whether the 
conservation of the grain number in the toppling rules is a 
necessary condition to obtain a critical state. It has been 
shown already that too much non-conservation leads to 
avalanches of characteristic sizes However, if grains 
are taken out of the system slowly, the system is found to 
be critical in some situations. A non-conservative version 
of the Abelian sandpile model with directional bias shows 
a mean field type critical behaviour |3^. Therefore, the 
detailed role of the conservation of the grain numbers 
during the topplings is still an open question. 

We acknowledge D. Dhar with thanks for a critical 
reading of the manuscript and for useful comments. 
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